Plant-derived compounds effectively inhibit the main protease of SARS-CoV-2: An in silico approach

The current coronavirus disease 2019 (COVID-19) pandemic, caused by the coronavirus 2 (SARS-CoV-2), involves severe acute respiratory syndrome and poses unprecedented challenges to global health. Structure-based drug design techniques have been developed targeting the main protease of the SARS-CoV-2, responsible for viral replication and transcription, to rapidly identify effective inhibitors and therapeutic targets. Herein, we constructed a phytochemical dataset of 1154 compounds using deep literature mining and explored their potential to bind with and inhibit the main protease of SARS-CoV-2. The three most effective phytochemicals Cosmosiine, Pelargonidin-3-O-glucoside, and Cleomiscosin A had binding energies of -8.4, -8.4, and -8.2 kcal/mol, respectively, in the docking analysis. These molecules could bind to Gln189, Glu166, Cys145, His41, and Met165 residues on the active site of the targeted protein, leading to specific inhibition. The pharmacological characteristics and toxicity of these compounds, examined using absorption, distribution, metabolism, excretion, and toxicity (ADMET) analyses, revealed no carcinogenicity or toxicity. Furthermore, the complexes were simulated with molecular dynamics for 100 ns to calculate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), and hydrogen profiles from the simulation trajectories. Our analysis validated the rigidity of the docked protein-ligand. Taken together, our computational study findings might help develop potential drugs to combat the main protease of the SARS-CoV-2 and help alleviate the severity of the pandemic.

Introduction plausible anti-SARS-CoV-2 effects according to an in silico prognosis [29]. Likewise, emodin, hesperidin, and chrysin have shown an S protein inhibiting efficacy that is analogous to chloroquine in tandem with hydroxychloroquine [30]. Phytochemicals such as lignan, quercetin, kaempferol, N-cis-feruloyltyramine, sugiol, jasmonic acid, putaminoxin D, 5,7-dimethoxyflavanone-4 0 -O-β-d-glucopyranoside, bonducellpin D, and caesalmin B have displayed anti-SARS-CoV-2 effect against Mpro target [31][32][33]. Moreover, karuquinone B, Lonicerae Japonicae Flos, and castanospermine, against S protein, in conjunction with cryptotanshinone, moupinamide, quercetin, coumaroyltyramine, kaempferol, tanshinone IIa, and N-cisferuloyltyramine, against PLpro protein, have revealed inhibitory effects on SARS-CoV-2 [32,34,35]. Recently, a flavonol named quercetin, an FDA-approved compound included in antiallergy and antioxidant medicines, has shown a promising inhibitory effect against SARS--CoV-2. As natural compounds and phytochemicals are efficiently processed at lower costs and have an impact on every stage of the interaction between the host and the virus, the quest for the most efficacious anti-SARS-CoV-2 therapeutic phytochemical is a worthwhile investigation to be performed both in silico and in vivo. Therefore, natural products and several phytochemicals, solely or in combination with traditional therapies, can be utilized to preclude and treat SARS-CoV-2 disease [29].
Computational drug detection schemes that include virtual screening techniques and molecular dynamics simulation proceedings are reliable high-throughput methods to identify plausible antiviral phytochemicals among a diverse array of repurposed phytochemical candidates. We hypothesized that in conjunction with molecular docking, a centralized virtual screening scheme can estimate the binding energy of compatible molecular interactions between the objective protein substrate and the potential ligand library. Moreover, computational methods identifying possible binding modes and binding site residues located in conserved motifs are plausible and scientifically reputable techniques to recognize repurposing phytochemicals. Such methods using drug-repurposing phytochemicals with antiviral activities consume less time than laborious conventional screening methods.
In low-income, economically backward, and underdeveloped countries, 17.8% of the people have received only one dose of the COVID-19 vaccine (https://ourworldindata.org/covid-vacc inations?country). The lack of vaccines and effective antiviral components and their limited efficacy necessitate the development of phytochemical-based eco-friendly phytopharmaceuticals against viral diseases that inhibit viral replication and penetration, while having manageable side effects and cost-effectiveness [36,37]. Therefore, in this study, we aimed to identify effective inhibitors and therapeutic targets for blocking the function of SARS-CoV-2. We retrieved 1154 phytochemicals by literature review and docked them against Mpro using a molecular docking technique that calculates the binding affinities and modes between the target substrate (e.g. proteins) and numerous ligands, such as phytochemicals, in the shortest time possible.

Protein preparation
The Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB) Protein Data Bank was used to obtain the crystal structure of Mpro in SARS-CoV-2 (PDB-ID: 6LU7; resolution 2.16) [38]. The latest version of the Discovery Studio [39] and the PyMol software [40] package were used to deplete and capacitate the extracted protein structure through the computational process. First, we eliminated all the heteroatoms, water molecules, and inhibitors in Mpro. Then, a computational tool named SWISS PDB Viewer software [41] package was applied for further purification and energy depreciation in the presence of the Groningen Molecular Simulation (GROMOS) 43B1 force field. Alongside, the optimized crystal conformation of Mpro was scrutinized using the same software to recognize subsistent characteristics that include side-chain geometry, lacking hydrogen, and erroneous bond order.

Ligand preparation
About 1154 phytochemicals were identified from various medicinal plants after an extensive and thorough literature survey (S1-S3 Tables). The three-dimensional structures of the procured phytochemicals that would be efficient ligand molecules were retrieved from the Pub-Chem database [42]. For optimization and energy minimization of the ligand structure within about 2000 minimization steps, mmff94 force field [43] along with the steepest gradient algorithm was utilized.
Molecular docking study. Molecular docking through PyRx Virtual Screening Tools (Autodock Vina, v.1.2.0.) [44] was performed to determine the binding affinity between the phytochemicals and Mpro and to elucidate the binding pose, exhibiting every single feasible orientation and conformation for any specific ligands at the binding site of the Mpro and the phytochemicals. A suitable gradient method and a universal force field (UFF) were used to achieve structural optimization and energy minimization. The numerical value of the sum of the number of minimization steps programmed was about 2000. As compatible ligands and substrates are required to undergo docking in AutoDock Vina, Mpro was chosen as the macromolecule substrate, and phytochemicals were selected as the ligands. After precisely preparing the ligand and the substrate, a grid box in Autodock was imposed to identify the substratebinding pocket representing the active site of the main protease. At the docking site of Mpro, the assigned grid box annotated the dimensions of X: 50.3334 Y: 67.2744 Z: 59.2586 Å, centered on X: 26.299 Y: 12.6039 Z: 58.9455 Å. As the ligand-flexible and protein-fixed docking was run, an anchored protease was identified, where every bond of the ligands was rotatable. After the completion of molecular docking, the highest docking energy signified the preeminent conformation. Based on the highest binding affinities and non-bonded interactions, the three top-most potential phytochemicals were marked for further investigation. Next, we used BIOVIA Discovery Studio and PyMol to explore the transcendent or docked conformations with the maximum docking energy by utilizing and assessing the visualized non-bonded interactions between the selected phytochemicals and the Mpro protein. Furthermore, the molecular docking approach with the formerly formulated grid was executed to dock the cocrystallized N3 ligand of 6LU7 by utilizing the 'PyRx' tool.

Absorption, distribution, metabolism, excretion, and toxicity (ADMET) analysis
Indispensable online-based servers, such as the ADMET structure-activity relationship (admetSAR) [45], SwissADME [46], and pKCSM [47] databases, were used to analyze and evaluate the physicochemical descriptors in conjunction with the pharmacokinetic properties. The above-mentioned web servers require Canonical Simplified Molecular-Input Line-Entry System (SMILES) retrieved from the PubChem database that predicts medicinal chemistry compatibility of the selected plausible antiviral phytochemicals.

Biological activities of the drug candidates
The online cheminformatics server "Molinspiration" (https://www.molinspiration.com/) was used to evaluate the biological activity of the screened phytochemicals. SMILES of these phytochemicals retrieved from the PubChem database were supplied as inputs to assess several biological activity parameters. This program applies sophisticated Bayesian statistics to assess a training set that consists of active structure and compares with the inactive molecules. Based on the investigation, a fragment-based model was created, and the bioactivity of each substructure fragment was estimated, adding up to the total activity contributions of the molecules' fragments.

Molecular dynamics simulation
The structural stability and constancy of the docked complex were precisely evaluated through the molecular dynamics simulation study. For molecular dynamics simulation, another scientific artificial reality application (YASARA) dynamics program [48] which used Assistant Model Building with Energy Refinement 14 (AMBER14) forcefield [49] was utilized. A cubic simulation cell was generated and extended for 20 Å at every part of the simulated proteinligand complex. The docked complexes were initially cleaned and optimized with the hydrogen bond network orientations. The TIP3P water solvation model was used at 0.997 g/L-1, 25˚C, and 1 atm, and energy minimization was performed using the steepest gradient algorithms with simulated annealing methods [50]. The physiological condition of the simulation cell was neutralized by the addition of 0.9% NaCl at pH 7.4 and 25˚C. Long-range electrostatic interactions were calculated using the Particle Mesh Ewalds algorithm with a cut-off radius of 8 Å [51]. The time step of the simulation was set as 1.25 fs. The simulation trajectories were saved after every 100 ps, and the simulation was extended up to 100 ns periods. The trajectories were used to calculate the root mean square deviation (RMSD), root mean square fluctuation (RMSF), radius of gyration (Rg), solvent-accessible surface area (SASA), and hydrogen bonds [52][53][54][55][56][57][58][59][60][61].
Our investigation revealed that cosmosiine interacted with Mpro through three conventional hydrogen bonds at positions Thr190, Gln189, and Glu166, in addition to three pi-alkyl bonds formations at positions Met165, Met49, and Cys145 (Fig 2 and Table 1). On the other hand, pelargonidin-3-O-glucoside formed six conventional hydrogen bonds at Glu166, His172, Leu141, Asn142, His163, and Thr190, one pi-pi T-shaped bond at His41, and one pialkyl bond at MET49. Besides, a multiplex of cleomiscosin A and Mpro was equalized by four conventional hydrogen bonds at Cys145, Glu166, Gly143, and Ser144, one alkyl bond at Met165, in addition to a pi-alkyl bond at Pro168. After docking the co-crystallized N3 ligand with the prepared grid, the binding affinity was -7 kcal/mol, which was comparatively lower than the previously encountered binding affinity such as -8.4, -8.4, and -8.2 kcal/mol. Although the N3 co-crystallized ligand resided in the binding pocket of 6LU7, the lower binding affinity regarding N3 interaction reveals the preferentiality of 6LU7, substantiating virtual screening in tandem with docking analyses.

ADMET
Pharmacokinetics and toxicity, alongside other drug-like properties, of the screened top-most phytochemicals, were assessed for evaluating the potency and safety (Table 2). Moreover, several operative properties of these top phytochemicals, including carcinogenicity, hepatotoxicity, p-glycoprotein inhibition, cytochrome P450 (CYP) inhibition, and central nervous system (CNS) permeability, were evaluated. Here, CNS permeability refers to a compound's ability to pass through the selectively semipermeable blood-brain barrier [47,62]. It has been shown that the CNS permeability value should be greater than -2 to penetrate the central nervous system [63]. Neither toxicity nor carcinogenicity is subsistent within the selected phytochemicals. The molecular weight (MW) of the screened molecules, cosmosiine, pelargonidin-3-O-glucoside, and cleomiscosin A, were 432.4 g/mol, 433.4 g/mol, and 386.4 g/mol, respectively. These values

PLOS ONE
Inhibition of SARS-CoV-2 main protease by plant compounds indicate that these compounds adhere to the range of MW for the Lipinski rule of five [64]. Six donors and 10 acceptors of hydrogen bond in cosmosiine, seven donors and nine acceptors of hydrogen bond in pelargonidin-3-O-glucoside, and two donors and eight acceptors of hydrogen bond in cleomiscosin A were observed. Besides, each of these hit phytochemicals showed neither hepatotoxicity nor acute oral toxicity and followed the Lipinski rule of five.

Biological activities of the drug candidates
Several observations required meticulous evaluation to assess the biological activities of the top three plausible phytochemicals with significant antiviral repurposing therapeutic potential (Table 3). Cosmosiine showed the most preeminent GPCR ligand activity when compared with pelargonidin-3-O-glucoside and cleomiscosin A. In the case of ion channel inhibitor activity, pelargonidin-3-O-glucoside exhibited a comparatively lower value than cosmosiine and cleomiscosin A. As a kinase inhibitor, cleomiscosin A manifested an increased activity than pelargonidin-3-O-glucoside, whereas cosmosiine showed the most predominant activity.
Cosmosiine also behaved as a more potent nuclear receptor-ligand compared to the other two phytochemicals. As a protease inhibitor, cosmosiine revealed the highest activity over both pelargonidin-3-O-glucoside and cleomiscosin A. All the screened phytochemicals exhibited promising enzyme inhibitor activities with cosmosiine having the highest activity value, and pelargonidin-3-O-glucoside and cleomiscosin showing almost similar activity values.

The molecular dynamics simulation study
The root mean square deviations of the simulations complexes were investigated to observe the variations and stability of the complexes . Fig 3(a) demonstrates how the top three docked complex, cosmosiine, pelargonidin-3-0-glucoside, cleomiscosin A showed an initial upward trend in RMSD. Cleomiscosin A maintained stability from the 20 ns time frame, and a lower degree of deviation in RMSD was observed for the entire simulation period. Pelargonidin-3-0-glucoside maintained a trend similar to that of cleomiscosin A for the maximum simulation periods. Pelargonidin-3-0-glucoside surpassed cleomiscosin A in RMSD in the 65-75-ns period, which might have ascended due to the higher instability of this complex. However, it quickly began to show a lower RMSD profile by maintaining a similar trend like cleomiscosin A, which also displayed a similar RMSD pattern like these two complexes in the entire simulation period. It had a slightly higher deviation at the last 85-90-ns time, although it did not significantly deviate. The solvent-accessible surface area (SASA) of the simulated complexes was measured to determine how the protein volume changed, with a greater SASA profile indicating increased protein surface area and a lower SASA profile indicating truncation of the protein complexes. Fig 3(b) demonstrates that the cleomiscosin A complex had a higher SASA profile than the other two complexes in the whole simulation trajectories, which correlated with the complexes' extensions. The expansion of the protein volume in SASA might have occurred due to the breakage of internal bonding. The other two complexes, with cosmosiine and pelargonidin-3-0-glucoside, had lower but constant SASA profiles in the whole simulation trajectories.
The radius of gyration of the simulation complexes is related to their labile nature, where a higher Rg profile indicates a more mobile nature of the complexes . Fig 3(c) shows that the three top complexes had a steady Rg profile in the entire trajectory and did not deviate Note: Bioactivity score > 0 (biologically active); −5.0 < Bioactivity score < 0 (moderately active); Bioactivity score < 0 (biologically inactive). https://doi.org/10.1371/journal.pone.0273341.t003

PLOS ONE
Inhibition of SARS-CoV-2 main protease by plant compounds significantly. The lower degree of fluctuations in the Rg descriptors for all three complexes correlated with the constant nature of the docked complexes. We also analyzed the hydrogen bond from the trajectories as it had a vital role in determining the stability of the docked complexes. The top three screened complexes had regular

PLOS ONE
Inhibition of SARS-CoV-2 main protease by plant compounds hydrogen bonding patterns in the simulating environment, with cleomiscosin A having the most number of hydrogen bonds [Fig 3(d)].
The docked complexes were again analyzed after 100 ns simulation and explored for their binding interactions to evaluate their binding rigidity. The cosmosiine and the Mpro complexes were stabilized by five hydrogen bonds at Gln192, Gln189, Arg188, Thr26, and Thr190 positions, and three hydrophobic bonds (two pi-alkyl and one pi-sulfur bond) at Met49, Met165, and Cys145 positions ( Table 4). The pelargonidin-3-O-glucoside and the Mpro complexes exhibited four hydrogen bonds at Glu166, His163, Ser144, and Leu141 after the 100 ns simulation time. Cleomiscosin A complexes had three hydrogen bonds at Gly143, Cys145, Asn142, and two alkyl bonds at Met165 and His41 positions.

Discussion
According to our analysis and screening via molecular docking and dynamics investigations, three potent phytochemicals were selected, which showed higher binding affinity for Mpro at their active sites than the other compounds, which is essential for the targeted inhibition of Mpro. In the case of cosmosiine, after the 0 ns simulation time, we observed one interaction in Domain 1 (Met49), three interactions in Domain 2 (Cys145, Met165, and Glu166), and two interactions in the extended loop region (Gln189 and Thr190) connecting Domain 2 and Domain 3 (Fig 2 and Table 1). After the 100-ns simulation time, the cosmosiine and Mpro complex exhibited two interactions at Domain 1 (Thr26 and Met49), two interactions at Domain 2 (Met165 and Cys145), and four interactions at the more extended loop region (Gln192, Gln189, Arg188, and Thr190) which connects Domain 2 and 3. Cosmosiine can inhibit the cytotoxicity and the dysregulation associated with ACE2, IL1α, and TGFβ expressions induced by the recombinant spike protein of SARS-CoV-2 [65]. It has also been found to have several biological activities, such as induction of apoptosis, autophagy, and cell cycle arrest [66]. Previous studies have also suggested that cosmosiine inhibits cell migration in cervical cancer [67], insulin receptor phosphorylation, adiponectin secretion [68], inflammation, and oxidative stress [69]. Moreover, cosmosiine has potential anti-inflammatory and antibacterial properties [69].
According to a previous study by Abha et al., 2011, cleomiscosin A had a more significant binding capability with toll-like receptors (TLR-4), cluster of differentiation molecules (CDs), and inducible nitric oxide synthase (iNOS) protein [75]. There have been no studies on its interaction with the Mpro of SARS-CoV-2. In our study, cleomiscosin A displayed six interactions at Domain 2 (Cys145, Glu166, Gly143, Met165, Pro168, and Ser144) after the 0 ns simulation time (Fig 2 and Table 1). The stabilization of cleomiscosin A and Mpro complexes manifested one interaction at Domain 1 (HIS41), in conjunction with four interactions at Domain 2 (Gly143, Cys145, Asn142, and Met165) after the 100 ns of simulation time (Table 4). Cleomiscosin A also exhibits anti-inflammatory, analgesic, and antipyretic properties [76], as well as antihepatotoxic and anti-inflammatory properties [77,78].
During the entire simulation period, diverse common interactions concerning the simulation sets for cosmosiine were identified, including Met49, Met165, Cys145, Gln189, and Thr190, which indicates the binding stability of the complex throughout the simulation time. Constant interactions were observed at Domain 1 (Met49), Domain 2 (Met165, and Cys145), and the linker loop region (Gln189 and Thr190). Our findings showed that cosmosiine predominantly interacts with the active site of the Mpro enzyme through residues including Met165, Cys145, Gln189, and Thr190, among which Cys145 is highly conserved among all coronaviruses [14]. Furthermore, pelargonidin-3-O-glucoside displayed distinct substantial interactions in the simulation epoch (during the 0 ns to 100 ns simulation time), including Glu166, His163, and Leu141. Precisely, the static interactions were remarkable at Domain 2 (Glu166, His163, and Leu141). However, pelargonidin-3-O-glucoside interacts with Mpro at the leading site through the residues Glu166, His163, and Leu141. Several stable interactions within the entire simulation time were identified, including Gly143, Cys145, and Met165 in the case of cleomiscosin A; stable interactions were observed at Domain 2 (Gly143, Cys145, and Met165). Prominently, cleomiscosin A interacts with Mpro enzyme's active site through residues including Gly143, Cys145, and Met165. These common interactions throughout the simulation period reveal the binding rigidity or the binding stability of the complex.
The cosmosiine and pelargonidin-3-O-glucoside compounds exihibit beta-D-glucopyranosyl moeity where the pelargonidin-3-O-glucoside functions as a plant metabolite and cosmosiine has role in non-steroidial drug, metabolite and antibacterial agents. Both ligand molecules are beta glucosidage and functions in metabolisams. The cleomiscosin A is a heterocyclic compound which exihibit antiinflammatory activity. This compound is beta lactone and organic heterocyclic compounds.
Molecular dynamics simulation of the docked complexes also revealed that the three compounds, cosmosiine, pelargonidin-3-O-glucoside, and cleomiscosin A, had a stable profile in several simulated trajectories, including RMSD, RMSF, SASA, Rg, and hydrogen bond studies. We also acquired images from various simulated trajectories, including at 25, 50, 75, and 100 ns, to investigate for changes in the binding sites or the binding rigidity at the binding pockets (Figs 4-6). In the simulated images, the docked positions of the three complexes were stable. The findings of the molecular docking and dynamics simulation studies suggest that these three compounds potentially inhibit the function of the SARS-targeted CoV-2 Mpro. However, these findings must be verified in a wet-lab setting.

PLOS ONE
Inhibition of SARS-CoV-2 main protease by plant compounds

Conclusion
Using a computational method, we identified effective inhibitors of the SARS-CoV-2 main protease. The phytochemical library was tested against the main protease using molecular docking to identify the most potent lead compounds. In addition, three of the top ligand molecules in the library, namely, cosmosiine, pelargonidin-3-O-glucoside, and cleomiscosin A, were found to bind to the active site of the targeted protein. The binding orientation and the stiffness of the docked protein-ligand complexes were also determined using molecular dynamics simulation. Simulation descriptors, such as RMSD, RMSF, SASA, and Rg, and hydrogen bond descriptors aided in analyzing the inflexible character of the complexes in atomistic environments. The toxicity and carcinogenicity of the top compounds were thoroughly investigated using multiple computational tools, and a possibility of no adverse or unfavorable effects was observed. Because this study relies solely on computational mining, these findings require further in silico cross docking and in vitro investigations, including enzymatic assays, to confirm the computational findings.
Supporting information S1